%% Load the data - Select the File: 'ViscosityProfiles\0pt5PEO'
close all; clear all; clc;
load('.\ViscosityProfiles\0pt5PEO')

vis = []; vis = cell2mat(mu');
xlim_max = 1000; x = []; x = xPos + (500 - median(xPos));

FitViscSlope = []; first = 4; last = 8;
for i = 1:length(mu)
    clear Mu_Exp p
    Mu_Exp = flip(mu{1,i}(first:last)); p = polyfit(xPos(first:last),Mu_Exp,1);
    FitViscSlope(i,1) = p(1);FitViscSlope(i,2) = p(2);
    FitVisc = polyval(FitViscSlope(i,:),xPos(first:last)); % Create line of fitted data points
end
Grad_Eta = mean(FitViscSlope(4:10,1));

st = 1; fin = 10;
LinVisProfile = figure('color',[1 1 1])
    col = 256*length(st:fin);
    cmapdef = colormap(jet(256*length(st:fin)-255)); %Define Colormap
    col_int = [256 linspace(256,col-256,length(st:fin)-1)+256]-255;
    colmap = cmapdef(col_int, :);
    for i = 1:length(st:fin)
        plot(x,vis(i,:),'color',colmap(i,:));
        hold on;
    end
    hold off;

    ax = gca;
    ax.FontSize = 7;
    xlim([0 xlim_max]);ylim([0 4.5]); xticks([0 250 500 750]); yticks([0,1,2,3,4,5]);
    set(gca, 'FontName', 'Arial');
    xlabel('Position x (\mum)','FontSize',8);ylabel('Viscosity \eta (cP)','FontSize',8);
    % Color Bar
    cbar_label = (st:fin)*interval - interval;
    colormap(jet);
        h = colorbar('color',[0 0 0]);
        h.Ticks = linspace(0,1,length(st:fin));      
        h.LineWidth = 1
        h.TickLabels = num2cell(cbar_label);
        h.FontSize = 7;
        h.FontName = 'Arial'; h.Location = 'eastoutside';
        title(h,'Time (min)')
% Plot MSD of each bin
Seq = 2;
Seq_Save = num2str((Seq-1)*10);
min = 1; maxTime = 25; t = 1:maxTime; t_s = t./fps;

col = 256*(CellSize(2)+2); cmapdef = colormap(copper(256*(CellSize(2)+2)-255)); %Define Colormap
close all;
col_int = [256 linspace(256,col-256,(CellSize(2)+2)-1)+256]-255;
cmap = cmapdef(col_int, :); leg = {};

PlotMSD_Adjy_SinTime = figure('color',[1 1 1]);
    p = zeros(1,CellSize(2));
    for j = 1:CellSize(2) % bins to evaluate
        leg{j} = num2str(x(j));
        p(j) = plot(t_s(min:2:maxTime),MSD_y_adj{1,Seq}(j,min:2:maxTime),...
            'color',cmap(j+1,:),'LineStyle', 'none','marker','o','LineWidth',1,'MarkerFaceColor',[1 1 1]);
        hold on;
    end
    fit = []; time_fit = []; fit_y = [];
    for j = 1:CellSize(2) % bins to evaluate
        fit(j,:) = polyfit(t_s(min:maxTime),MSD_y_adj{1,Seq}(j,min:maxTime),1);
        t_s_new = [t_s 7];
        y1 = polyval(fit(j,:),t_s_new);
        f = plot(t_s_new,y1,'-','color',cmap(j+1,:),'LineWidth',1);
        uistack(f,'bottom');
        time_fit = t_s_new; fit_y(:,j) = y1;
    end
    
    x = round(x,0); legLabel = [0]; legLabel = [legLabel x 1000];
    colormap(copper);
        h = colorbar('color',[0 0 0]);
        h.Ticks = linspace(0,1,CellSize(2)+2);    
        h.LineWidth = 1;
        h.TickLabels = num2cell(legLabel);
        h.FontSize = 7;
        h.FontName = 'Arial'; h.Location = 'eastoutside';
        %title(h,'Bin Center')
    
    ax = gca; ax.FontSize = 7;
    xlim([0 6.5]);ylim([0 11.5]);
    xticks([0 2 4 6]); yticks([0 5 10]);
    set(gca, 'FontName', 'Arial');
    xlabel('Time (s)','FontSize',8);ylabel('MSD (\mum^2)','FontSize',8);
 
% Crop Figure
    alw = 1; % Border line width
    afs = 7; % Font Size
    markersize = 2;
    xlabel(''); ylabel('');
    mypapersize = [3.9 2.9]; %[width height]
    set(gca, 'FontName', 'Arial')
    set(gca,'LineWidth',alw);
    set(gca,'XColor',[0 0 0]);set(gca,'YColor',[0 0 0]);
    set(gcf, 'PaperUnits', 'inches', 'papersize', mypapersize, 'PaperPosition', [0 0 mypapersize]);
    set(PlotMSD_Adjy_SinTime,'Units','inches','Position',[0 0 mypapersize]);

% Save Data
dir_Excel = '.\Source data\ExtendedData\';
% Experiment
tabA = 1;
writematrix(t_s(min:maxTime)',[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','A4')
writematrix(x,[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','C2')

writematrix(MSD_y_adj{1,Seq}(1,min:maxTime)',[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','C4') % 77 um
writematrix(MSD_y_adj{1,Seq}(2,min:maxTime)',[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','D4') % 161 um
writematrix(MSD_y_adj{1,Seq}(3,min:maxTime)',[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','E4') % 246 um
writematrix(MSD_y_adj{1,Seq}(4,min:maxTime)',[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','F4') % 331 um
writematrix(MSD_y_adj{1,Seq}(5,min:maxTime)',[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','G4') % 415 um
writematrix(MSD_y_adj{1,Seq}(6,min:maxTime)',[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','H4') % 500 um
writematrix(MSD_y_adj{1,Seq}(7,min:maxTime)',[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','I4') % 585 um
writematrix(MSD_y_adj{1,Seq}(8,min:maxTime)',[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','J4') % 669 um
writematrix(MSD_y_adj{1,Seq}(9,min:maxTime)',[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','K4') % 754 um
writematrix(MSD_y_adj{1,Seq}(10,min:maxTime)',[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','L4') % 839 um
writematrix(MSD_y_adj{1,Seq}(11,min:maxTime)',[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','M4') % 923 um
% Fit
writematrix(time_fit',[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','O4')
writematrix(x,[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','Q2')

writematrix(fit_y(:,1),[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','Q4') % 77 um
writematrix(fit_y(:,2),[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','R4') % 161 um
writematrix(fit_y(:,3),[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','S4') % 246 um
writematrix(fit_y(:,4),[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','T4') % 331 um
writematrix(fit_y(:,5),[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','U4') % 415 um
writematrix(fit_y(:,6),[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','V4') % 500 um
writematrix(fit_y(:,7),[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','W4') % 585 um
writematrix(fit_y(:,8),[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','X4') % 669 um
writematrix(fit_y(:,9),[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','Y4') % 754 um
writematrix(fit_y(:,10),[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','Z4') % 839 um
writematrix(fit_y(:,11),[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabA,'Range','AA4') % 923 um
% Save Figure
figdir = '.\Supplemental\ExtendedDataFigure5\'
print(PlotMSD_Adjy_SinTime, '-painters','-dpdf', '-r1200', [figdir,'\MSDvsTime_',Seq_Save,'Min.pdf']);

% Error Bars for Selected Viscosity Profile Plot
close all;clc
x = xPos + (500 - median(xPos));
col = 256*(CellSize(2)+2);
cmapdef = colormap(copper(256*(CellSize(2)+2)-255)); % Define Colormap
close all;
col_int = [256 linspace(256,col-256,(CellSize(2)+2)-1)+256]-255;
cmap = cmapdef(col_int, :);
leg = {};

kT = 4.11*10^-21;   % Thermal Energy [J=kg*m^2/s^2]
D = 0.5;            % diameter [microns]
a = 1.0*10^-6*D/2;  % radius [meters]
Diffusion = {};eta = {}; Diffusion_EP = {}; eta_EP = {};

for j = 1:length(MSD_y_adj{1,Seq}(:,1))
    catMSDy = []; p = []; p_err = [];D = []; D_STD = [];
    catMSDy = 0.5*vertcat(MSD_y_adj{1,Seq}(j,:));
    [p,p_err] = LinearFitError(t_s(min:maxTime),(catMSDy(min:maxTime))); % p(1) = Diffusion Coefficient [um^2/s]
    D = p(1) * (1.0*10^-6).^2; D_STD = p_err(1) * (1.0*10^-6).^2; % Diffusion & Prop. of Error [m^2/s]
    eta{1,j} = (kT./(6.*pi.*a.*D))/1E-3 % [cP]
    eta_EP{1,j} = sqrt((-kT./(6.*pi.*a) * (1/D)^2 * D_STD)^2) / 1.0E-3 % [cP]
    
    Diffusion{1,j} = D; % Diffusion [m^2/s]
    Diffusion_EP{1,j} = D_STD; % Propagation of Error: Diffusion [m^2/s]
end

FillColor = colmap(Seq,:); LineColor = colmap(Seq,:);
LineWidth = 1;
alpha_errorbar = 0.20;
x = xPos + (500 - median(xPos));
y = cell2mat(eta);
error = cell2mat(eta_EP);
ViscProfile_fillErrorBar = figure('color',[1 1 1]);
        patch = fill([x,fliplr(x)],[y + error,fliplr(y - error)],FillColor);
        set(patch, 'edgecolor', 'none'); set(patch, 'FaceAlpha', alpha_errorbar);
        hold on;
        plot(x,y,'-','Color',LineColor)
        
    for i = 1:length(vis(Seq,:))
        plot(x(i),vis(Seq,i),'color',cmap(i+1,:),'color',cmap(i+1,:),'LineStyle', 'none','marker','o','LineWidth',1,'MarkerFaceColor',cmap(i+1,:));
    end
    hold off;
    x = round(x,0);
    legLabel = [0];legLabel = [legLabel x 1000];
    colormap(copper);
        h = colorbar('color',[0 0 0]);
        h.Ticks = linspace(0,1,CellSize(2)+2);     
        h.LineWidth = 1;
        h.TickLabels = num2cell(legLabel);
        h.FontSize = 7;
        h.FontName = 'Arial'; h.Location = 'eastoutside';
        %title(h,'Bin Center')
    
    ax = gca;
    ax.FontSize = 7;
    xlim([0 1000]);ylim([0 4.5]); xticks([0 250 500 750 1000]); yticks([0,1,2,3,4,5]);
    set(gca, 'FontName', 'Arial')
    xlabel(['Position x (',char(181),'m)'],'FontSize',8);ylabel('Viscosity \eta (cP)','FontSize',8);
    % Crop Figure
    alw = 1; % Border line width
    xlabel(''); ylabel('');
    mypapersize = [3.9 2.9]; %[width height]
    set(gca, 'FontName', 'Arial')
    set(gca,'LineWidth',alw);
    set(gca,'XColor',[0 0 0]);set(gca,'YColor',[0 0 0]);
    set(gcf, 'PaperUnits', 'inches', 'papersize', mypapersize, 'PaperPosition', [0 0 mypapersize]);
    set(ViscProfile_fillErrorBar,'Units','inches','Position',[0 0 mypapersize]);
% Save Data
dir_Excel = '.\Source data\ExtendedData\';
% Experiment
tabB = 2;
writematrix(x',[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabB,'Range','A3')
writematrix(cell2mat(eta)',[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabB,'Range','B3')
writematrix(error',[dir_Excel,'ExtDataFigure5.xlsx'],'Sheet',tabB,'Range','C3')

% Save Figure
figdir = '.\Supplemental\ExtendedDataFigure5\'
print(ViscProfile_fillErrorBar, '-painters','-dpdf', '-r1200', [figdir,'\ViscProfileErrorBar_',Seq_Save,'Min.pdf']);
%% FUNCTION: 
function [p,p_err] = LinearFitError(x,y)
% This function calculates the error of the fit to a linear line: Y = BX + A
    N = length(x);
    delta = N*sum(x.^2)-sum(x)^2;

    A = (sum(x.^2)*sum(y) - sum(x)*sum(x.*y))/delta; % intercept 
    B = (N*sum(x.*y) - sum(x)*sum(y))/delta;         % slope
    p = [B A];                                       % mirror format of polyfit in matlab

    sigma_y = sqrt(1/(N-2) * sum((y-A-B*x).^2));

    sigma_A = sigma_y*sqrt(sum(x.^2)/delta); % intercept error
    sigma_B = sigma_y*sqrt(N/delta);         % slope error
    p_err = [sigma_B sigma_A];
end